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The recently developed feedback trap can be used to create arbitrary virtual potentials, to explore 
the dynamics of small particles or large molecules in complex situations. Experimentally, feedback 
traps introduce several finite time scales: there is a delay between the measurement of a particle's 
position and the feedback response; the feedback response is applied for a finite update time; and a 
finite camera exposure integrates motion. We show how to incorporate such timing effects into the 
description of particle motion. For the test case of a virtual quadratic potential, we give the first 
accurate description of particle dynamics, calculating the power spectrum and variance of fluctuations 
as a function of feedback gain, testing against simulations. We show that for small feedback gains, the 
motion approximates that of a particle in an ordinary harmonic potential. Moreover, if the potential 
is varied in time, for example by varying its stiffness, the work that is calculated approximates that 
done in an ordinary changing potential. The quality of the approximation is set by the ratio of the 
update time of the feedback loop to the relaxation time of motion in the virtual potential. 

PACS numbers: 05.40.-a, 87.19.1r, 87.15.Vv 



I. INTRODUCTION 

In 2005, Cohen and Moerner introduced the Anti- 
Brownian ELectrokinetic (ABEL) trap, a new experimen- 
tal technique for studying long-time dynamical properties 
of small particles and molecules. One key advantage over 
other trapping techniques such as optical or magnetic 
tweezers is its ability to trap molecules and sub-micron 
particles directly rather than via the micron-sized particles 
of the former techniques. The ABEL trap uses feedback to 
counteract the random thermal fluctuations that perturb 
the motion of small objects in a finite-temperature fluid 
PQ, with electrokinetic forces being one way among many 
to apply restoring forces; it is thus perhaps more simply 
termed a feedback trap. The basic idea is to observe the 
position of an object, compare its estimated position with 
a desired position, and, as rapidly as possible, apply a 
corrective force to move the particle towards the desired 
position. To trap small objects, the observations must be 
rapid, as the particle diffuses "out of control" during the 
time t s between corrections. Indeed, the lower limit to 
the size of particle that can be trapped depends directly 
on t s . Recently, Fields and Cohen, by responding to every 
detected photon, were able to trap for several seconds a 
single fluorescent dye molecule diffusing in water [2]. 

Cohen has also shown that a feedback trap can do 
more than just trap a particle: it is also possible to place 
the particle in an arbitrary "virtual" potential [3L The 
protocol to approximate motion in a potential U(x) works 
as follows: Let the estimated position of the particle at 
time step n, in units of a sampling period t s , be x n . (We 
distinguish between the observed position x n and the true 
position x n .) Then, at time step n, we apply the force 
F„ = —d x U(x n ), held constant over the interval t s . In 
[3], experimental evidence is given that the probability 
distribution of observed positions x, obeys the Boltzmann 
distribution p(x) oc exp[— U{x)/kBT]. 

In the above scheme, the potential imposed is a virtual 



one that is imposed by the rules of the feedback loop. 
Take away the feedback loop, and there is only a particle 
diffusing in a fluid. But then it is fair to ask, In what sense 
is the motion of a discrete closed-loop feedback system 
equivalent to a "real" potential? After all, the dynamics 
is discretized at a time scale t s . Reasoning by analogy to 
computer simulations of the Langevin equation, we expect 
that as t s — > 0, the closed- loop feedback system will be 
equivalent to the desired continuous dynamical system. 
But experiments are always done at finite t s , and the in- 
formation is acted upon only after a delay, t^. In addition, 
the position is typically measured (with random error) 
by integrating a camera over a finite time t c that can be 
comparable to t s . Under such conditions, is the motion 
truly equivalent to the desired potential? Are there correc- 
tions to the "naive potential" ? If so, are such corrections 
important in a given application? Likewise, can one use 
virtual potentials for thermodynamic calculations such as 
the work done by a changing potential? 



We will find that the answers to these questions are 
mathematically boring — the situation resembles that of a 
continuous system with discrete observations — but physi- 
cally exciting: we can now study the motion of a particle 
in an arbitrary potential or force field and learn about 
both dynamic and thermodynamic quantities. 



Indeed, this article was motivated by attempts to use a 
feedback trap to explore Landauer's Principle that erasing 
information in a memory element requires a finite amount 
of work [H [5] . After making preliminary measurements 
reporting a qualitative observation of the effect [BJ, we 
realized that there were systematic deviations from calcu- 
lations based on a continuous potential that needed to be 
understood before a quantitative study could be made. 
The calculations reported below address these concerns. 
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II. DYNAMICS OF A PARTICLE IN A 
QUADRATIC VIRTUAL POTENTIAL 

In this section, we explore the dynamics of a feedback 
trap with an imposed quadratic virtual potential. The 
goal will be to calculate the power spectrum and variance 
of the particle's position fluctuations. The calculation is 
complicated by the presence of three short time scales, 
which are comparable but not in general equal: 

1. The output is updated after an observation with a 
delay time td . 

2. The observation is made via a camera exposure of 
duration t c . 

3. The update is applied for a time t s , which is also 
the periodicity of camera exposures (one exposure 
starts every t s , which may be longer than t c ). 

As we will see, previous work has not described all the 
consequences of these experimental complications. Below, 
we will see how to incorporate them by solving a series of 
increasingly complicated problems: 

1. The diffusion dynamics of a Brownian particle with 
observations every t s that are integrated over t c , 
with a response that is delayed by a time td = t s ; 

2. feedback trapping in a virtual quadratic potential, 
with noiseless, instantaneous position measurements, 
whose results are available with no delay; 

3. feedback trapping adding a camera exposure t c and 
a delay t d — t s ; 

4. feedback trapping in the general case where td ^ t s . 

A. Free diffusion 
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FIG. 1. Timing diagram for feedback trap with instantaneous 
position measurements. The feedback update interval is t s , 
and there is a delay td = t s between the measurement time 
and the output of the next feedback update. The true position 
of particles indicated at bottom, together with the 

associated thermal noise £ n and force F n . The force F n is 
computed from the observed position x n available at the time 
F n is started, nt a . At top are indicated the observed position 
x n and the amount of associated thermal noise £ n in each 
observation. 



fluctuation-dissipation theorem [HE], (t)) = and 
(fM (£)£(")(£')) = 2DS(t - t'), where the diffusion con- 
stant D — ksT/j and 5(t — t') is the Dirac delta function. 
Integrating the continuous equation of motion over the 
time interval [nt s , (n + l)t s ) then gives 

X n +1 = X n + £„ , f« = / C (t,) (*) dt . (2) 

Jnt, 

where represents the displacement due to thermal 
forces integrated over time interval n. Again, £ n are Gaus- 
sian random variables with mean and with (£ m £n) — 
2Dt s 5 mn , where 8 mn is the Kronecker delta function. 

We next consider the distinction between the actual 
position x n and the observed position x n that becomes 
available at the same time. Taking into account that the 
observed position is based on an exposure of duration t c 
whose midpoint is delayed by t d = t s , we have 



We begin by considering the one-dimensional free diffu- 
sion of a Brownian particle that is observed via camera 
exposures that measure the average position over the 
exposure time t c . There is one exposure every t s inter- 
val, and the mid-point of the exposure time is delayed 
by one time step, td = t s . We distinguish between the 
actual position of a particle at time nt s , denoted by x n 
and the corresponding observation of that position (as 
deduced from the camera exposure), denoted x n . The 
above statements are illustrated in the timing diagram 
shown in Fig. [TJ 

We start with the equation of motion of the Brownian 
particle, assuming complete overdamping and no applied 
forces: 

x = -^(t) = ^\t), (1) 
7 

where ^ F \t) is the fluctuating thermal force and ^ v \t) 
gives the corresponding velocity fluctuations. From the 



x n +i = — / x(t)dt = x n - f n 0) + £„ , 

tc Jnt a — ^t c 

rnt s 

^(t)dt, 

Jnt,-\t c 

e„ = - / dt / ^\t')dt'. (3) 

tc Jnt s ~±t c Jnta — ^tc 

In Eq. ([3]), x n — £ n (0) is the position of the particle at the 
beginning of the exposure. The term £ n represents the 
thermal noise as averaged by the camera exposure and 
gives the average displacement from the position at the 
beginning of the exposure to its end. The negative ^i -* 
is needed because we define the measured position S n +i 
relative to the position at the midpoint of the camera 
exposure and not at the beginning. (This convention will 
be convenient for what follows.) 

Next, we consider the statistics of the measured particle 
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displacement 



= £n-l 
= Cn-l 



with similar definitions for all other quantities [x = Z(x n ), 
etc.). Neglecting initial conditions, we then have 
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Because the thermal noise terms in Eq. ^ are all Gaus- 
sian with mean zero, we immediately have that (Ax n ) = 0, 
which simply says that there is no bias to the displace- 
ments of a random walker. In the last line of Eq. Q, 



Sn-1 — Sn-1 Sn + 



t(0) 
' Sn-1 



nt s — 



(n-l)t. 



(5) 



which is just the thermal noise over one sample period, 
f s , starting and ending at the beginning of the camera 
exposure. By contrast, £n-i is the thermal noise over one 
period starting and ending at the midpoint of the camera 
exposure. 

The mean-square displacement during t s is then 

((Ax n ) 2 ) = ((£„-i) 2 )+ 2 ((£„) 2 )- 2<C-i4Vi) • (6) 

In Eq. we omit the terms 2{£ f n _ 1 £ n ) = -2(£„_i £„} = 
(no overlap). We can evaluate 
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Multiplying Eq. ( 10b ) by z — 1 and substituting Eq. ( 10a 
gives 

(z - l)zx = £-(z- 1)^ (0) + (z- l)f 

= ? + {z-i)i, (ii) 

where we make the same redefinition of the thermal noise 
term that we did in Eq. |5]), with (,'(z) the Z-transform 
of the thermal noise £„. The power spectrum (modulus 
of discrete-time Fourier transform) is then 

|z(z-l)| 2 <|*| 2 ) = 

2[<ie'i 2 ) + k-ii 2 (ieT) + ^-i)<ro+ c - c -] . (12) 

where c.c. denotes the complex conjugate of the last terms 
and where the factor of 2 results from considering only 
positive frequencies. Evaluating the thermal noise expres- 
sions and substituting z = e ltJJta , we have 

\z(z - 1)| 2 (\x\ 2 ) = 2[2Dt s + 2(1 - cos ujt s )(lDt c ) 
+ (e iut > - l)(Dt c ) + c.c] 
= 4Dt s - \Dt c {l - cosut s ) . (13) 
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Z (v \t')^(t")) dt" Solving for (|x| 2 ) as a function of the angular frequency 



we have 
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where we have shifted the domain of integration for all 
the integrals by (n — l)t s — |i c , for clarity. A similar 
calculation by Cohen [9] gives ((£n) 2 ) = \Dt c . Putting 
these results together, we have 



((Ax n ) 2 ) =2Dt s + 2(%Dt, 
= 2D(t s 



■ 2 



2{Dt c ) 



(8) 



B. Harmonic potential, with "perfect" 
measurements 

We begin our study of virtual potentials in the simplest 
case, where the measurements are "perfect": that is, they 
are instantaneous (t c = 0), without delay (td = 0), and 
free of observation noise. We maintain an update time 
of t s . Although no experiment is so simple, the results 
already illustrate some of the key differences between 
motion in a true and in a virtual potential. Under these 
assumptions, the equations of motion Eq. ([2| become 



which gives the finite-exposure-time correction to the 
usual mean-square displacement [9TfTT]. 

Turning now to the power spectrum of the position 
measurements, we take the Z-transform (or equivalcntly, 
calculate the generating function) of Eqs. ([2| and Let 
us define the ^-transform of the sequence x n to be 



kra+l 



(15) 



We first note that Eq. ( 15 ) is identical to the Euler algo- 



rithm for integrating stochastic differential equations |12j . 
In Eq. (15), we can view the a term as deriving from a 



Z[x n ] = x{z) = ^2 



(9) 



n=0 



"proportional feedback" law, F n = —kx n , with a — t s k/j. 
Although the k in the feedback law superficially resembles 
the force constant k^ in Hooke's law, F(t) = —k^x(t), 
it is not quite the same quantity. To see this point qualita- 
tively, we note that the force in Hooke's law continuously 
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changes as x(t) changes, whereas the force in the feedback 
system F n is constant over the interval t s . As long as 
d«l, the difference between the two situations is not 
great. But larger values of a lead to different dynamics. 

To illustrate this point, we calculate the steady-state 
variance. Assuming (x 2 n+1 ) — = (x 2 ), squaring 



Eq. (15 1, and noting that (x n ^ n ) — 0, we have 



(* 2 > = (l-a)V) + (a, 



(16) 



where (&) = (£ 2 ) = 2Dt s 
we then have 



With 1 - (1 - a) 2 = a(2 - , 



(x 2 ) = 



2DU 



a(2 - a) 



(17) 



For a <C 1, the variance is Dt s /a, which is the value 
expected by a "naive" application of the equipartition 
principle: substituting a — t s k/^ and D = UbT/^ gives 
(x 2 ) = ksT/k. However, the variance at finite a is always 
larger and diverges at a* = 2, beyond which the motion is 
unstable. Physically, the extra variance comes from over- 
correcting perturbations. At a = 2, the motion oscillates 
from one side of equilibrium to the other. 

Repeating the argument for unit delay, td = t s , we have 



Xn+l X Ti 



(18) 



For this case, squaring and averaging Eq. ( 18 ) leads to 

-a 2 {x 2 ) 



2Dt, 



2a(x X-i) 



evaluated by multiplying Eq. (18| by x n and averaging, 
giving (xx—i) = (x 2 ) — a(«_i] 



The last term can be 
by x n and a 1 
The result is 



= 2Dt s 



1 



a(l - a){2 + a) 



(19) 



which is unstable for a > 1 and again goes to the equipar- 
tition result Dt s /a when a — > 0. The two expressions for 
(x 2 ), along with corresponding simulations, are plotted as 
a function of feedback gain a in Fig. |2j We note that the 
reduction in critical gain a* from 2 to 1 reflects the perfor- 
mance deterioration caused by the delay. Longer delays 
further decrease a* . For example, a similar calculation 
gives a* = \{\fb - 1) 0.62 for t d = 2t s . 



C. Harmonic potential, with camera exposure 

We next include a finite camera exposure t c and also 
observation noise %„. The equations of motion gener- 
alize in two ways. First, we define the force F n to act 
between time nt s and (n + l)t s . As Fig. [I] shows, the 
force changes midway through the camera exposure time. 
Integrating the equation of motion for the continuous 
dynamics forwards and backwards in time from t n and 
neglecting stochastic forces, we find, 



x{t) 



c n Fj^—xitji t~j t n —\ <c t <c t n 

&n ~t~ —F n (t — tfi) t n <C t <C t n -\-i 



(20) 
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FIG. 2. (Color online.) Position variance for a virtual harmonic 
potential in a feedback trap, for t s = 0.01 s and D = 2 /im 2 /s. 
Simulations based on Eqs. (15| and ( |18| l are for the same 
parameters, for a total simulation time of r = 1000 s. The 
solid round markers correspond to the case o f td = 0, and the 
accompanying theory curve is given by Eq. (171. The hollow 
round markers correspond to the case td = ts, and the accom- 
panying theory curve is given by Eq. | ]19[ ). The dashed line 
shows the "naive" variance calculated from the equipartition 
theorem (x 2 ) = (feT/fe) = Dt s /a. The vertical dashed gray 
line indicates unstable motion, which occurs beyond a* = 2 
for td = and a* = 1 for td — t s . 



Inserting the expression for x(t) in Eq. (20) into the 
average for x n+ i in Eq. p]), we have 



l-n+l 



-F n -i(t n — t) 



dt 



x n T F n (t ^n) 
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(21) 



Equation (21 ) implies that the measured position is biased 
by a difference in forces in intervals n and n + 1. In 
particular, there will be no bias if the forces are constant, 
F n _i — F n - In that case, there is a uniform drift; the 
midpoint of the exposure is at x n ; and the particle moves, 
on average, an equal amount before and after the midpoint 
(in time) of the exposure. 

Taking into account the unit delay td — t s leads to 
coupled, linear difference equations for x n and x n - 

Xn+l — X n OLXn T £n 

X n+ i = X n + Oi{x n -l - X n ) - £, n 0) + £n + Xn , (22) 

where a = is the ratio of the sampling time to the 

trap relaxation time, a' = a-^- measures the effects of 
the finite camera exposure time t c , and where we have 
included stochastic terms in the x equation. In particular, 
Xn accounts for observation noise due to photon shot noise, 
microscope resolution, and other effects. For simplicity, 
we assume it to be Gaussian, with (x n ) = and (x 2 ) = x 2 - 
The distinction between the state variable (the position 
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x n ) and its measurement x n in Eq. (|22|) is emphasized in which implies 



discussions of control theory [131 H4j. 

To find the power spectrum, we generalize slightly the 
analysis of Sec. |II A| Taking the Z-transform, we have 

(z-l)x = -ax + £, (23) 



(z - l)zx = -ax - a' ^- ) x + £ + (z - 1)(£ + X ) 



[z 2 - (1 - a')z + {a - 2a') + a'z^x = £' + («- 1)(£ + x ) . (24) 
Solving for x then leads to an expression for the power spectrum: 

c + (z-m+ x ) 



[z 2 - (1 - a')z + (a - 2a') + a 1 z 



'r-H 



(25a) 



= 2[2Dt s + 2(x 2 -|^ e )(l~cos^)] 
Ml/ | e 2 Iwts _ (i _ a /) e io>t s + ( a -2a') + a'e-^ t °\ 2 ' 1 J 



To our knowledge, Eq. ( |25[ ) has not been previously 
derived. Previous versions [Sj [TS] neglect the a' terms 
in the denominator. Physically, those terms are present 
because the averaged position x n and not x n is used 
in the feedback loop. In general, the denominator in 
transfer functions such as Eq. (25a I reflects the structure 
of feedback loops 



loop. 



In Fig. [3j we illustrate typical power spectra for low 
(a = 0.1) and high (a = 0.9) values of the feedback 
gain. Markers denote simulations using parameter values 
for D, t c , td, t s , and x that are typical of experimental 
systems. The a = 0.1 case approximates the Lorentzian 
spectrum expected of a continuous system, although even 
in that case, the discrete observations ensure significant 
deviations from the continuous spectrum for frequencies 
comparable to the Nyquist frequency (50 Hz here). The 
dashed line shows the corresponding Lorentzian spectrum, 
which, in our notation, is given by 



ADti 



AD 



+ {iut s y (i/t?) + < 



(26) 



On the other hand, the high-gain curve shows a signifi- 
cant resonant peak. Physically, the peak corresponds to 
an overshoot. Unlike a continuous quadratic potential, 
the force applied in a feedback trap is constant over the 
interval t s . For large feedback gains, the particle can 
overshoot the set point, leading to decaying damped os- 
cillations and a resonance in the power spectrum. The 
corresponding Lorentzian spectrum (dashed line) agrees 
with our calculation at low frequencies but then excludes 
the power associated with the resonance. The instability 
threshold (a* « 1.14) is slightly larger than for t c = 0, 
where a* = 1. The camera exposure acts as a low-pass fil- 
ter of the position measurement, stabilizing the feedback 
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FIG. 3. (Color online.) Power spectra for a virtual harmonic 
potential in a feedback trap, for t s — 0.01 s, (td/t s ) = 1, 
(tc/ts) = 0.95, D = 2 fim 2 /s, x = 0.018 ^m, and a = 0.1 and 
0.9. Simulations are for the same parameters, for simulation 
time t = 400 s. The power spectra (solid lines) are plotted 
from Eq. ( |25b[ ) and are not fit to the simu latio ns. The corre- 
sponding Lorentzian approximations, Eq. (261, are shown as 
dashed lines. 



In Fig. |4j we illustrate the simulated and predicted 
variance of position measurements as a function of the 
feedback gain a. Although the variance could be calcu- 
lated using strategies similar to those used above for the 
t c = case, it is simpler to integrate the power spectrum 
in Eq. (25b I from to the Nyquist frequency (2/t s ). For 
a <C 1, the "naive variance" predicted from equipartition 
is a good approximation to the exact value. At a = 0.1, 
the delay increases the variance by about 15%. 
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FIG. 4. (Color online.) Position variance for a virtual har- 
monic potential in a feedback trap, for t a = 0.01 s, (td/t 3 ) = 1, 
(t c /t s ) = 0.95, D = 2 ^m 2 /s, x = 0.018 ^m. Simulations 
(hollow round markers) are for the same parameters, for sim- 
ulation time r = 400 s. The predicte d valu es (solid black 
curve) are evaluated by integrating Eq. ( 25b I over frequency 
and are not fit to the simulations. The dashed line shows the 
"naive" variance calculated from the equipartition theorem 
(a; 2 ) = (ksT/k) = Dt s /a. The vertical gray line indicates 
unstable motion for a > a* ~ 1.14. 



D. Harmonic potential, general delay 

With zero delay, the discrete equation of motion for 
a particle in a feedback trap with harmonic potential is 



a first-order difference equation, Eq. (15 1. A delay of 
(td/t s ) = 1, in turn, leads to a second-order difference 



equation, Eq. (18). Since an nth-order difference equation 
is the discrete analog of an nth-order differential equation, 
does a fractional delay, (td/t s ) ^ integer, lead to an 
analog of a fractional derivative? Such derivatives have 
an analytic structure that is qualitatively different from 
that of integer derivatives [16] . In fact, the answer is 
much less exotic, and the results for fractional delay are 
qualitatively similar to those for integer delay. 

To see this, we will assume that td ~ t s . In particular, 
we will assume, as illustrated in the modified timing dia- 
gram, Fig.[5j that the measured position x n +i is computed 
from a camera exposure that straddles the contributions 
from two forces, F n -i and F n . To generalize to completely 
arbitrary time delays requires a somewhat awkward no- 



tation that separates the number of integer periods of t s 
contained in t d and its remaining fractional part and also 
distinguishes between the straddling case treated here 
and the non-straddling case (similar to and simpler than 
the one treated here). 
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FIG. 5. Timing diagram for feedback trap with instantaneous 
position measurements. The timing is similar to Fig. [I] except 
that the delay time td is not equal to the sampling time t s . 

The first step is to repeat the derivation of the bias in 
x n +i- The difference from Eq. (21 ) is that the integration 
limits now are [t n +i~t d -\t c , tn)< and [t„, t n+1 —td+~t c ). 
We now find 



x n +i = x n + !r~ (a+F„ - a-F n -i) , (27) 
07 



where 



a± = 



1± 



2(t s - t d ) 



(28) 



Note that a± — 1 when td — t s , and Eq. ( 27 1 reduces to 



Eq. (21 1. The apparent divergence at t c — > is an artifact 
of our assumption that the exposure straddles two forces, 
which is not true in this limit. 

The coupled equations for x n and x n become 



Xn+1 — X n dX n <vn 

Xn+1 — X n ~\~ Oi ((2_X n _i Q+Xn 



ei 0) + 1, 



(29) 



Taking the Z-transform then gives the transfer function 
and power spectrum: 



^(^X^*) (30a) 



{z 2 — (1 — a'a + )z + [a — a'{a + + a_)] + a'a^z 1 } 



/|-,2X = 2[2^ s + 2( X 2 -|^ e )(l-cos^)] 

^ | e 2^t 3 _ (1 _ a i a+ )e iu}t ° + [a - a'(a+ + a_)] + a'a_e- lL ^ | 2 ' 1 ; 



As promised, the power spectrum and related results III. VIRTUAL POTENTIALS AND 

are just slightly more complicated than for integer delay. THERMODYNAMICS 
Indeed, it is generally true that fractional delays merely 

change coefficients in a discrete dynamical system [13]. Can virtua i potentials be used for thermodynamic cal- 

culations? In this section, we investigate the accuracy 
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of "naive" calculations of the work done by a changing 
potential, following ideas of stochastic thermodynamics 
[T7HT5] . As an example calculation, we calculate the mean 
work required to vary the stiffness of a virtual harmonic 
potential in a finite time. We will find that estimates of 
work agree, to 0(a), with those of a true potential. 



Harmonic potential with varying force constant 

To explore the work done by a virtual potential, we con- 
sider a time-dependent potential U(x,t). We start with 
the case of a quadratic virtual potential with a feedback 
gain a n that is increased in constant steps from a« at 
t = to a f at t = t = Nt s . Recall that a„ corresponds 
to a force F n = —k n x n , but k n only approximates the 
force constant of a harmonic potential. 

To begin, we recall the calculation of the average work 
done in the continuous case, where the true force constant 
fc( c )(i) varies from fc| to fey . We further assume that 
the variation is done slowly enough (r — > oo) that, at 
each moment, the system is in local equilibrium. Then, 



W 



dU 



dt 



- I k^x 2 ,H . 



Assuming that (x 2 )(t) = fc B T/[fc( c )(i)], we have 



(W) 
k B T 



fcW , 1 

—r-rdt = - 
fc( c ) 2 



dk^ 



= 1 ln ~ f 



(31) 



(32) 



For fixed time t, we define the partition function Z(t) = 
fT oo exp[-U(x,t)/k B T]dx ~ fc^ 1 / 2 . Then, in terms of 
the free energy J- = — k B T hi Z, we have simply 



(W) = AT , 



(33) 



as expected for an adiabatic protocol. The calculation in 
Eq. ( 32 ) can be generalized to finite r [20] . 



Consider next the equivalent calculation for a feedback 
system with td — t c — \ — (no delay, no camera 
exposure, no observation noise). Discretizing Eq. (31) 
gives 



k B T ~ 2 ^ k B T [Xn> s ~ 2D ^ 



O > (34) 



where d = (a.f — <Xi)/r is constant. Rewriting Eq. (17) as 

1 1 



x z = Dt„ 



a 2 — a 



(35) 



we have 



(W) _ t s 



N 



n— v 

where a n = on + ant s is the "staircase" feedback gain. 



k B T 2 ^ \a„ 2 — a r , 

n— v 



In the Appendix, we show that 



lim t„ 



T/t. 

n=0 



= In 



b + at 
b + cu. 



(37) 



Applying this result to Eq. (36) with 6 = and —2, we 
have 



(mi 

k B T 2 



a. 



a f 



(38) 



We can quickly generalize to the case of unit delay. 
Rewriting Eq. ( 19 ) as 



DU 



1 

o 



1/3 4/3 



a 



(39) 



leads to 



m 

k R T 



In 



a f 



1/3 



4/3 



(40) 

Equations |38|) and (40) agree with simulations. 

Figure [6] illustrates the work calculations and shows the 
approach to steady-state behavior for finite cycle time r. 
In Fig.[6ja), we see the approach to the value (horizontal 
black line) calculated from Eq. ( 40 1 . The solid curve is a 
fit, for large r (> 1), to the extra (dissipated) work 
that exceeds the asymptotic value calculated in Eq. (40). 
The fit is to the form Wd ~ t _1 , as suggested by Sekimoto 

In Fig. pi b) , we see that the distribution of work has 
a pronounced positive skew for finite transition times 
r but approaches a Gaussian distribution as r — > oo 
(as expected in general for continuous systems [21]). A 
positive skew in the work distribution has been seen before 
in experimental measurements based on optical tweezer 
traps |19j , but a quantitative theory is so far not available. 
It is interesting that the skew predicted here, which should 
be accessible experimentally, is much greater than that 
observed in [T^] . 

In Fig.^c), we see that the skew decreases as r -1 / 2 for 
large r. To reduce the effect of rare outliers, we use a ro- 
bust measure of skew, [(W75 — W 50 ) — (W 50 — W 2 5))/W 5 o 7 
where W50 is the median value of work and W25 and VF75 
represent the first and third quartile values [22]. Skew- 
ness thus serves to measure disequilibrium. An alternative 
measure — requiring more data to calculate with equiva- 
lent precision — is based on the relative entropy between 
forward and backward processes |23j . The origin of the 
scaling of the skew measure deserves more investigation. 

We note that the values for the asymptotic work differ 
from the "naive" calculation of 4 In by O(a) for a -C 1. 
Since the work depends on (x 2 ) and since this quantity 
also differs by 0(a) from the continuous value for small 
a, the observation is expected. The main point is that 
there are no secular terms whose error grows with r. 
Rather, no matter how long the measurement, the errors 
are controlled by the discretization size a. 



(a) 
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FIG. 6. (Color online.) Work required to vary a harmonic 
potential, with a varying between a; = 0.05 and Q/ = 0.4 
for te. = t s = 0.01s and D = 2 /jm 2 /s. (a) Average work as 
a function of transition time r. Solid line is a fit to a t" 1 
correction for r > 1 s. (b) Histogram of work measurements 
for t = 1, 8, and 100 s. (c) Robust skewness vs. t. 



We also note that the work we calculate is greater than 
that for the corresponding continuous case. In a recent 
paper, Sivak et al. have pointed out that the integra- 
tion scheme itself injects extra "shadow" work into the 
system [24]. In principle, integration schemes that are 
more sophisticated than Euler integration could be used 
to remove the effects of shadow work. 

Finally, we have explored numerically the effects of 
a camera exposure t c and observation noise \ an d find 
that they have only a small numerical effect (data not 
shown) . Such a result is expected from the results of our 
calculations of power spectra and variance in Section [TTJ 
where we saw that t c and \ made only minor shifts in 
those quantities. 



IV. DISCUSSION AND CONCLUSIONS 

In this article, we have shown that the complications 
brought on by finite time scales in virtual potentials can 
be accounted for and do not lead to radical differences 
between virtual and actual potentials, as long as we are 



careful. We find that calculations of dynamics — such as 
variance, power spectrum, and work — lead to answers 
that agree with continuum calculations to 0(a), where a 
is the ratio of the feedback update time to the relaxation 
time of the potential. The situation is then similar to that 
of ordinary experiments that also discretize continuous 
signals. The main practical difference is that experiments 
with a physical potential can usually sample at very small 
values of a (e.g., a — 3 x 10~ 4 in [25 J, whereas feedback 
experiments, which must act in real time on each observa- 
tion, often use higher values (e.g., a = 0.1 in J]). As we 
saw in Section [II B| minimizing the delay helps to increase 
the value of the instability gain a* and, by extension, to 
increase the range of acceptable feedback gains a. 

Although feedback traps use relatively large values of 
a, they can impose more general and better-controlled 
motion (via virtual potentials) on a particle than can be 
done with physical potentials. For example, in explor- 
ing the Landauer principle with the feedback trap, we 
have been able to implement the scheme proposed by 
Dillenschneider and Lutz [25], which is based on a spe- 
cific parametrization of a double-well potential where the 
difference in well depth and barrier heights are separately 
controlled. Although such manipulations can be approx- 
imately done by placing two optical traps side by side 
and adding a sideways flow [5], more general and better- 
controlled variations are possible in the feedback trap, 
since the form of the imposed potential is now arbitrary, 
as long as the feedback gain is small enough. 

One assumption in the work presented above is that 
the electrokinetic force in an ABEL trap responds instan- 
taneously to the command signal. This is usually a good 
approximation: both electrophoretic and electroosmotic 
forces depend on a realignment of the double-layer, which 
is typically much less than a micron in scale. Typical 
response times are then of order 10 /is , which is much 
faster than the ms time scale of feedback traps. In related 
work, we have used piezoelectric translation stages to 
track the motion of small particles rather than trap them. 
The formalism that is required is similar to that explored 
here, but one must then account for the finite relaxation 
time, of order ms, of the translation stages [27l [28] , 

We have also shown that calculations of work done 
by the virtual potential are relatively straightforward, 
although many theoretical questions deserve more ex- 
ploration. However, the situation is less clear regarding 
the calculation of heat dissipation into a thermal bath. 
A naive discretization analogous to that used for work, 
Eq. ( 34 ) gives the correct result when applied to harmonic 



potentials but has spurious secular terms that grow lin- 
early with the protocol time for anharmonic potentials. 
However, we observe similar effects when using the Euler 
discretization to calculate heat dissipation for ordinary 
Langevin equations with anharmonic potentials (but no 
feedback). Indeed, how to calculate the heat properly 
from sampled data — with or without feedback — remains 
an issue of current discussion jTHl EH EH] • 

One final question is whether experiments can impose 
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conditions and forces that are well-enough controlled to 
apply the theory described in this work. That is, we have 
assumed that the particle moves only because of thermal 
fluctuations and because of deliberately applied forces. 
Knowing their effects requires a calibration of both elec- 
trical mobility (the charge on a particle) and diffusion 
constant (the size of a particle). It also assumes that the 
electric field is known and there are no mechanical drifts. 
In preliminary experimental work, we have established 
that it is possible to accurately describe motions at fre- 
quencies above 1 Hz, but at lower frequencies, mechanical 
drifts from temperature variations become important. 



Proof. Intuitively, the expression is simply J^ f jq^- 
More formally, with Aa = a/ — a; and N — r/t s , we 
have 



N 

E 



N ^b- 

n— 



f 



Aan 
N 



N 

E ~(i>+ . , A 
n=0 Aa 



Aa 



n 
Aa 



(42) 
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APPENDIX 



We prove Eq. ( 37 ) , which is used in the calculation of 



The digamma function ip(x) ~ In a; for x —> oo. Taking 
N — » oo, we have 



In „ A ° , « In ' Aq 



(b+aj)N / I (b+og) 

Aa I \ Aa 

/ A, 
= In I 



b + a t 
b + aj 



average work. = m j __T L j _ (43) 

Lemma. 

^tf^f^iYy 1 =ln ( b +^l). (4i) 

r^oo \ T J ' b + a n \ b + 04 J 
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